The aim of exploration seismology is to obtain an image of the subsurface by probing it with seismic waves at various locations. These waves are generally generated by using airguns in marine seismics, and vibroseis or dynamite in land seismics. The waves propagate downwards through the subsurface, and are reflected at interfaces between geological layers or refracted within layers. Parts of these waves subsequently propagate upwards to the surface, where they are detected and recorded.
In exploration seismology, although the time coordinate is regularly sampled, spatial coordinates are typically irregularly sampled due to the presence of obstacles in land and strong currents in marine. But even receivers placed for example within a marine survey cable or streamer may not be always equidistant. Hence, the inline sampling in direction of the cable can be quite irregular.
The interpolation on regular grid points or regularization of seismic signals or data is very important, particularly for uses in time-lapse survey matching, multiple suppression and imaging. If the irregular nature of the sampling grid is ignored or handled poorly, notable errors are introduced, the severity of which may be further amplified at later stages of the seismic processing chain.
The problem of signal reconstruction from uniformly spaced data has been investigated in depth.
The Whittaker-Kotel'nikov-Shannon sampling theorem states that any signal f(x) can be reconstructed from its uniformly spaced samples, if the sampling interval is less than half the period of the highest spectral component in that signal. Thus if f(x) is bandlimited to the wavenumber σ/2, which is known as the Nyquist wavenumber, then the sampling theorem provides the following formula to interpolate any function value from uniformly spaced values k(m/σ):
                              f          ⁡                      (            x            )                          =                              ∑                          m              =                              -                ∞                                      ∞                    ⁢                                    f              ⁡                              (                                  m                  /                  σ                                )                                      ⁢            sin            ⁢                                                  ⁢                          c              ⁡                              (                                                      σ                    ⁢                                                                                  ⁢                    x                                    -                  m                                )                                                                        (        1        )            where sin c(x)=sin(πx)/(πx). Thus, when the sampling rate is sufficient and there is no aliasing, the sampling theorem provides a way to reconstruct the signal “exactly” from its samples uniformly spaced samples. To satisfy requirements of the sampling theorem, the signal should be sampled at a rate greater than twice the Nyquist rate, i.e., σ.
Moreover, if the signal is space limited, i.e. the signal has only negligible amount of energy outside an interval [0, (L−1)/σ], the infinite summation in Whittaker-Kotel'nikov-Shannon sampling theorem can be replaced with a finite summation:
                              f          ⁡                      (            x            )                          =                              ∑                          m              =              0                                      L              -              1                                ⁢                                    f              ⁡                              (                                  m                  /                  σ                                )                                      ⁢            sin            ⁢                                                  ⁢                                          c                ⁡                                  (                                                            σ                      ⁢                                                                                          ⁢                      x                                        -                    m                                    )                                            .                                                          (        2        )            
A feature of the sampling theorem is that, the infinite summation given in (1) is uniformly convergent in x. This allows a user to differentiate the expression of eq. (1) term-by-term and derive the reconstruction theorems for derivatives of f(x). The general formula for rth derivative is:
                                          f                          (              r              )                                ⁡                      (            x            )                          =                              ∑                          m              =                              -                ∞                                      ∞                    ⁢                                    f              ⁡                              (                                  m                  /                  σ                                )                                      ⁢                                          ⅆ                r                                            ⅆ                                  x                  r                                                      ⁢            sin            ⁢                                                  ⁢                                          c                ⁡                                  (                                                            σ                      ⁢                                                                                          ⁢                      x                                        -                    m                                    )                                            .                                                          (        3        )            
For the first derivative, r=1, this leads to the following interpolation formula:
                                          f            ′                    ⁡                      (            x            )                          =                              ∑                          m              =                              -                ∞                                      ∞                    ⁢                      σ            ⁢                                                  ⁢                          f              ⁡                              (                                  m                  /                  σ                                )                                      ⁢                                          {                                                                            cos                      ⁡                                              (                                                                              π                            .                                                    ⁡                                                      (                                                                                          σ                                ⁢                                                                                                                                  ⁢                                x                                                            -                              m                                                        )                                                                          )                                                                                                            σ                        ⁢                                                                                                  ⁢                        x                                            -                      m                                                        -                                                            sin                      ⁡                                              (                                                  π                          ⁡                                                      (                                                                                          σ                                ⁢                                                                                                                                  ⁢                                                                  x                                  .                                                                                            -                              m                                                        )                                                                          )                                                                                                            π                        ⁡                                                  (                                                                                    σ                              ⁢                                                                                                                          ⁢                              x                                                        -                            m                                                    )                                                                    2                                                                      }                            .                                                          (        4        )            
In applications where samples of the function and its higher order derivatives are available simultaneously, the required sampling rate for exact reconstruction can be reduced. For instance, in the simplest case where multichannel data consist of both data and the corresponding first order derivative (gradient) it is possible to reconstruct the continuous signal from samples recorded at Nyquist rate (rather than twice the Nyquist rate) by using the following:
                              f          ⁡                      (            x            )                          =                              ∑                          m              =                              -                ∞                                      ∞                    ⁢                                    .                              {                                                                                                                              f                          ⁢                                                      (                                                          2                              ⁢                                                              m                                /                                σ                                                                                      )                                                                          +                                                                                                                                                                          (                                                      x                            -                                                          2                              ⁢                                                              m                                /                                σ                                                                                                              )                                                ⁢                                                                              f                            ′                                                    ⁡                                                      (                                                          2                              ⁢                                                              m                                /                                σ                                                                                      )                                                                                                                                              }                                      ⁢            sin            ⁢                                                  ⁢                                                            c                  2                                ⁡                                  (                                                            σ                      ⁢                                                                                          ⁢                                              x                        /                        2                                                              -                    m                                    )                                            .                                                          (        5        )            
The relation (5) is known as multi-channel sampling theorem.
In other words, to reconstruct the continuous signal it is sufficient to know either signal samples f(m/σ) at twice the Nyquist rate or signal and gradient samples f(2 m/σ) and f′(2 m/σ), respectively, at the Nyquist rate. Also, it is interesting to note that in the multichannel sampling theorem the sin c function is squared. There exists a corresponding version of the multichannel sampling theorem for cases when samples of both signal and its Hilbert transform are available.
By definition, a signal bandlimited to σ/2 does not have any spectral components above the wavenumber σ/2. Therefore filtering such a signal with an ideal low-pass filter with a cut off wavenumber of σ/2 does not change its spectral content. Hence the following identity holds:
                                                                        f                ⁡                                  (                  x                  )                                            =                                                                    f                    ⁡                                          (                      x                      )                                                        *                                ⁢                σ                ⁢                                                                  ⁢                sin                ⁢                                                                  ⁢                                  c                  ⁡                                      (                                          σ                      ⁢                                                                                          ⁢                      x                                        )                                                                                                                                          =                                  σ                  ⁢                                                            ∫                                              -                        ∞                                            ∞                                        ⁢                                                                  f                        ⁡                                                  (                          τ                          )                                                                    ⁢                      sin                      ⁢                                                                                          ⁢                                              c                        ⁡                                                  (                                                      σ                            ⁡                                                          (                                                              x                                -                                τ                                                            )                                                                                )                                                                    ⁢                                              ⅆ                        τ                                                                                                        ,                                                          (        6        )            where ‘*’ denotes the convolution operator and σ sin c(σx) is an ideal low-pass filter with a cut off wavenumber of σ/2. The Riemann sum approximation to this integral identity yields the following approximation, which is known as the sin c interpolator:
                                                        f              ⁡                              (                x                )                                      ≈                                          f                L                            ⁡                              (                x                )                                              =                      σ            ⁢                                          ∑                                  m                  =                  0                                                  L                  -                  1                                            ⁢                              Δ                ⁢                                                                  ⁢                                  x                  m                                ⁢                                  f                  ⁡                                      (                                          x                      m                                        )                                                  ⁢                sin                ⁢                                                                  ⁢                                  c                  ⁡                                      (                                          σ                      ⁡                                              (                                                  x                          -                                                      x                            m                                                                          )                                                              )                                                                                      ,                            (        7        )            where Δxm is the Jacobian weight, i.e., Δxm=xm+1−xm, and f(xm) is the value of the seismic data at irregular offset xm. It is important to note that, when Δxm=1/σ, the sin c interpolator is exact since
                              f          ⁡                      (            x            )                          =                              ∑                          m              =                              -                ∞                                      ∞                    ⁢                                    f              ⁡                              (                                  m                  /                  σ                                )                                      ⁢            sin            ⁢                                                  ⁢                          c              ⁡                              (                                  σ                  ⁡                                      (                                          x                      -                      τ                                        )                                                  )                                                                        (        8        )            by Whittaker-Kotel'nikov-Shannon theorem. On the other hand, when Δxm≠1/σ, the sin c interpolator provides only a crude approximation to the continuous signal. In this case, a better approach would be to invert (2) for the desired uniformly spaced signal values f(m/σ).
By using the matrix notation the equation (7) can be written as
                              h          =                                                    [                                                                                                    f                        ⁡                                                  (                                                      x                            1                                                    )                                                                                                                                                                        f                        ⁡                                                  (                                                      x                            2                                                    )                                                                                                                                                ⋮                                                                                                                          f                        ⁡                                                  (                                                      x                            L                                                    )                                                                                                                    ]                            ≈                                                [                                                                                                              s                          11                                                                                                                      s                          12                                                                                            ⋯                                                                                              s                                                      1                            ⁢                            L                                                                                                                                                                                        s                          21                                                                                                                      s                          22                                                                                            ⋯                                                                                              s                                                      2                            ⁢                            L                                                                                                                                                              ⋮                                                                    ⋮                                                                    ⋱                                                                    ⋮                                                                                                                                      s                                                      L                            ⁢                                                                                                                  ⁢                            1                                                                                                                                                s                                                      L                            ⁢                                                                                                                  ⁢                            2                                                                                                                      ⋯                                                                                              s                          LL                                                                                                      ]                                ·                                  [                                                                                                              f                          ⁡                                                      (                            0                            )                                                                                                                                                                                        f                          ⁡                                                      (                                                          1                              /                              σ                                                        )                                                                                                                                                              ⋮                                                                                                                                      f                          ⁡                                                      (                                                                                          (                                                                  L                                  -                                  1                                                                )                                                            /                              σ                                                        )                                                                                                                                ]                                                      =            Sg                          ,                            (        9        )            where σ/2 is the bandwidth of the signal f(x) and S is the sin c matrix with entries sij=sin c(σ(xi−j/σ)). If the matrix S is well conditioned than the seismic data at regular offsets can be computed by standard matrix inversion:g=S−1h.  (10)
Otherwise, a least squares minimum norm inversion should be preferred:g==(STW1S+W2)−1STW1h,  (11)where W1 usually chosen as a diagonal matrix whose mth diagonal entry is the Jacobian weight Δxm=xm+1−xm and W2 is usually chosen as a small multiple of identity matrix, i.e., W2=ε2I.
An interpolator in accordance with equation (11) has been described in: Yen J. L., 1956, On non-uniform sampling of bandwidth-limited signals, IRE Trans. Circuit Theory, CT-3, 251-257 (1956). It is therefore sometimes referred to as Yen's Interpolator of Type 1.
Many known interpolators used in seismics are variations of the Yen's interpolator.
The interpolators based on Yen's 1st theorem usually provide satisfactory results on non-aliased signals with little high-wavenumber content. However their performance degrades significantly when the interpolated signal has a substantial amount of high wavenumber spectral content.
Another shortcoming of the interpolators based on Yen's 1st theorem is that in order to solve equation (9), in general at least as many irregular sampling positions as regular sampling positions are required. Hence, if some seismic traces are dropped out, traces which reside at further locations must be used to solve the system of equations given by eq. (9). Usually this degrades the accuracy of the interpolated sample values.
Further, although Yen's 1st interpolator is exact for infinite length signals, it is an approximation when only a finite extent of the signal is available for interpolation.
The interpolation as discussed so far are non-adaptive techniques, where a block of irregular input data is interpolated to determine a block of regular output data. Once the matrix or in general the regularization kernel is computed, the same kernel can be used on all data recorded with the same sampling grid. Therefore these methods provide a fast solution when many recordings with the same irregularity are to be regularized. On the other hand, the existence of even a few aliased spectral components in the input data introduces large regularization errors. Hence, the non-adaptive techniques usually have no means to overcome aliasing without additional measurements.
One of the recently developed techniques to overcome the limitations of the non-iterative interpolation algorithm is the anti leakage Fourier transformation (ALFT), described for example by Sheng Xu et. al in: Geophysics Vol. 70, No. 4 (July-August 2005), P. V87-V95. The ALFT approach to seismic interpolation begins with the spatial discrete Fourier transformation of the seismic traces usually in the FX domain. To this purpose the non-uniform Fourier transformation is used, which is basically a Riemann summation approximation to the continuous-space Fourier transformation. In the second step of the ALFT, the most energetic spectral component of the non-uniform Fourier transformation is detected and the corresponding sinusoidal space domain signal is computed on the irregular grid by using inverse Fourier transformation. Assuming that the frequency and amplitude of the most energetic spectral component are accurately estimated, the estimated sinusoidal signal is then subtracted from the input data to prepare the input for the next iteration. Hence, the steps of ALFT can be summarized as:                1. Non-uniform Fourier transform        
      G    ⁡          (      k      )        =            ∑              m        =        1            L        ⁢          Δ      ⁢                          ⁢              x        m            ⁢              g        ⁡                  (                      x            m                    )                    ⁢              ⅇ                              -            j2                    ⁢                                          ⁢          π          ⁢                                          ⁢                      kx            m                                              2. Find dominant spectral component        
      k    *    =      arg    ⁢                  ⁢                  max        k            ⁢                                G          ⁡                      (            k            )                                                      3. First order signal estimate        
                    g        ^                    k        *              ⁡          (      y      )        =      Δ    ⁢                  ⁢    k    ⁢                  ⁢          G      ⁡              (                  k          *                )              ⁢          ⅇ              j        ⁢                                  ⁢        2        ⁢                                  ⁢        π        ⁢                                  ⁢                  k          *                ⁢        y                            4. Compute residual g(xm)=g(xm)−ĝk*(xm) and iterate        
However the known ALFT method is found to have a very slow convergence rate. Further the method is not immediately suited to handle multichannel data, for example data which include gradient measurements.
Given the problems of the existing interpolators it remains an object to find improved interpolators capable of interpolating data received by receivers at irregular locations to regular sampling locations.